#include <iostream>
#include <sys/time.h>
#include <arm_neon.h>
#include <pthread.h>
#include <semaphore.h>
#include <cstdlib>
using namespace std;

const int MAX_SIZE = 1024;
const int STEP = 64;
const int THREAD_COUNT = 3;

sem_t mainSemaphore;
sem_t workerStartSem[THREAD_COUNT];
sem_t workerEndSem[THREAD_COUNT];

float dataMatrix[MAX_SIZE][MAX_SIZE];

struct ThreadConfig {
    int threadIndex;
    int size;
};

void* WorkerRoutine(void* arg) {
    float32x4_t vectorA, vectorB, vectorC, vectorD;
    ThreadConfig* config = (ThreadConfig*)arg;
    int index = config->threadIndex;
    int size = config->size;

    for (int col = 0; col < size; ++col) {
        sem_wait(&workerStartSem[index]);

        for (int row = col + 1 + index; row < size; row += THREAD_COUNT) {
            vectorA = vld1q_dup_f32(dataMatrix[row] + col);

            int pos;
            for (pos = col + 1; pos + 3 < size; pos += 4) {
                vectorB = vld1q_f32(dataMatrix[col] + pos);
                vectorC = vld1q_f32(dataMatrix[row] + pos);
                vectorD = vmulq_f32(vectorA, vectorB);
                vectorC = vsubq_f32(vectorC, vectorD);
                vst1q_f32(dataMatrix[row] + pos, vectorC);
            }

            for (; pos < size; ++pos)
                dataMatrix[row][pos] -= dataMatrix[row][col] * dataMatrix[col][pos];

            dataMatrix[row][col] = 0.0f;
        }

        sem_post(&mainSemaphore);
        sem_wait(&workerEndSem[index]);
    }

    return nullptr;
}

void MatrixInitialization(int size) {
    for (int i = 0; i < size; ++i)
        for (int j = 0; j < size; ++j)
            dataMatrix[i][j] = static_cast<float>(rand()) / 10.0f;
}

void SerialMatrixDecomposition(int size) {
    for (int col = 0; col < size; ++col) {
        for (int pos = col + 1; pos < size; ++pos)
            dataMatrix[col][pos] /= dataMatrix[col][col];
        dataMatrix[col][col] = 1.0f;

        for (int row = col + 1; row < size; ++row) {
            for (int pos = col + 1; pos < size; ++pos)
                dataMatrix[row][pos] -= dataMatrix[row][col] * dataMatrix[col][pos];
            dataMatrix[row][col] = 0.0f;
        }
    }
}

void VectorizedMatrixDecomposition(int size) {
    float32x4_t vectorA, vectorB, vectorC, vectorD;
    for (int col = 0; col < size; ++col) {
        vectorA = vld1q_dup_f32(dataMatrix[col] + col);

        int pos;
        for (pos = col + 1; pos + 3 < size; pos += 4) {
            vectorB = vld1q_f32(dataMatrix[col] + pos);
            vectorC = vdivq_f32(vectorB, vectorA);
            vst1q_f32(dataMatrix[col] + pos, vectorC);
        }

        for (; pos < size; ++pos)
            dataMatrix[col][pos] /= dataMatrix[col][col];

        dataMatrix[col][col] = 1.0f;

        for (int row = col + 1; row < size; ++row) {
            vectorA = vld1q_dup_f32(dataMatrix[row] + col);

            int posIter;
            for (posIter = col + 1; posIter + 3 < size; posIter += 4) {
                vectorB = vld1q_f32(dataMatrix[col] + posIter);
                vectorC = vld1q_f32(dataMatrix[row] + posIter);
                vectorD = vmulq_f32(vectorA, vectorB);
                vectorC = vsubq_f32(vectorC, vectorD);
                vst1q_f32(dataMatrix[row] + posIter, vectorC);
            }

            for (; posIter < size; ++posIter)
                dataMatrix[row][posIter] -= dataMatrix[row][col] * dataMatrix[col][posIter];

            dataMatrix[row][col] = 0.0f;
        }
    }
}

void MultiThreadMatrixDecomposition(int size) {
    sem_init(&mainSemaphore, 0, 0);
    for (int i = 0; i < THREAD_COUNT; ++i) {
        sem_init(&workerStartSem[i], 0, 0);
        sem_init(&workerEndSem[i], 0, 0);
    }

    pthread_t threads[THREAD_COUNT];
    ThreadConfig configs[THREAD_COUNT];
    for (int i = 0; i < THREAD_COUNT; ++i) {
        configs[i].threadIndex = i;
        configs[i].size = size;
        pthread_create(&threads[i], nullptr, WorkerRoutine, &configs[i]);
    }

    for (int col = 0; col < size; ++col) {
        float32x4_t tempVector = vld1q_dup_f32(dataMatrix[col] + col);

        int pos;
        for (pos = col + 1; pos + 3 < size; pos += 4) {
            float32x4_t currentVector = vld1q_f32(dataMatrix[col] + pos);
            float32x4_t resultVector = vdivq_f32(currentVector, tempVector);
            vst1q_f32(dataMatrix[col] + pos, resultVector);
        }

        for (; pos < size; ++pos)
            dataMatrix[col][pos] /= dataMatrix[col][col];

        dataMatrix[col][col] = 1.0f;

        for (int i = 0; i < THREAD_COUNT; ++i)
            sem_post(&workerStartSem[i]);

        for (int i = 0; i < THREAD_COUNT; ++i)
            sem_wait(&mainSemaphore);

        for (int i = 0; i < THREAD_COUNT; ++i)
            sem_post(&workerEndSem[i]);
    }

    for (int i = 0; i < THREAD_COUNT; ++i)
        pthread_join(threads[i], nullptr);

    sem_destroy(&mainSemaphore);
    for (int i = 0; i < THREAD_COUNT; ++i) {
        sem_destroy(&workerStartSem[i]);
        sem_destroy(&workerEndSem[i]);
    }
}

void Compute() {
    timeval start, end;
    double duration;

    for (int size = STEP; size <= MAX_SIZE; size += STEP) {
        cout << "Matrix Size: " << size << endl;

        MatrixInitialization(size);
        gettimeofday(&start, nullptr);
        SerialMatrixDecomposition(size);
        gettimeofday(&end, nullptr);
        duration = (end.tv_sec - start.tv_sec) * 1000.0 + (end.tv_usec - start.tv_usec) / 1000.0;
        cout << "Serial Matrix Decomposition Time: " << duration << " ms" << endl;

        MatrixInitialization(size);
        gettimeofday(&start, nullptr);
        VectorizedMatrixDecomposition(size);
        gettimeofday(&end, nullptr);
        duration = (end.tv_sec - start.tv_sec) * 1000.0 + (end.tv_usec - start.tv_usec) / 1000.0;
        cout << "Vectorized Matrix Decomposition Time: " << duration << " ms" << endl;

        MatrixInitialization(size);
        gettimeofday(&start, nullptr);
        MultiThreadMatrixDecomposition(size);
        gettimeofday(&end, nullptr);
        duration = (end.tv_sec - start.tv_sec) * 1000.0 + (end.tv_usec - start.tv_usec) / 1000.0;
        cout << "Multi-Threaded Matrix Decomposition Time: " << duration << " ms" << endl;

        cout << "----------------------------------------" << endl;
    }
}

int main() {
    Compute();
    return 0;
}

/*
Matrix Size: 64
Serial Matrix Decomposition Time: 19.379 ms
Vectorized Matrix Decomposition Time: 12.873 ms
Multi-Threaded Matrix Decomposition Time: 83.2 ms
----------------------------------------
Matrix Size: 128
Serial Matrix Decomposition Time: 190.89 ms
Vectorized Matrix Decomposition Time: 106.753 ms
Multi-Threaded Matrix Decomposition Time: 196.756 ms
----------------------------------------
Matrix Size: 192
Serial Matrix Decomposition Time: 546.599 ms
Vectorized Matrix Decomposition Time: 433.834 ms
Multi-Threaded Matrix Decomposition Time: 475.654 ms
----------------------------------------
Matrix Size: 256
Serial Matrix Decomposition Time: 1186.75 ms
Vectorized Matrix Decomposition Time: 841.405 ms
Multi-Threaded Matrix Decomposition Time: 614.568 ms
----------------------------------------
Matrix Size: 320
Serial Matrix Decomposition Time: 2110.7 ms
Vectorized Matrix Decomposition Time: 1447.31 ms
Multi-Threaded Matrix Decomposition Time: 1406.37 ms
----------------------------------------
Matrix Size: 384
Serial Matrix Decomposition Time: 3748.35 ms
Vectorized Matrix Decomposition Time: 2860.52 ms
Multi-Threaded Matrix Decomposition Time: 2272.9 ms
----------------------------------------
Matrix Size: 448
Serial Matrix Decomposition Time: 6246.88 ms
Vectorized Matrix Decomposition Time: 4744.67 ms
Multi-Threaded Matrix Decomposition Time: 3471.55 ms
----------------------------------------
Matrix Size: 512
Serial Matrix Decomposition Time: 9508.35 ms
Vectorized Matrix Decomposition Time: 7164.05 ms
Multi-Threaded Matrix Decomposition Time: 5068.11 ms
----------------------------------------
Matrix Size: 576
Serial Matrix Decomposition Time: 13776.7 ms
Vectorized Matrix Decomposition Time: 10094.1 ms
Multi-Threaded Matrix Decomposition Time: 6939.54 ms
----------------------------------------
Matrix Size: 640
Serial Matrix Decomposition Time: 18620 ms
Vectorized Matrix Decomposition Time: 13553.1 ms
Multi-Threaded Matrix Decomposition Time: 9058 ms
----------------------------------------
Matrix Size: 704
Serial Matrix Decomposition Time: 24683.3 ms
Vectorized Matrix Decomposition Time: 18126 ms
Multi-Threaded Matrix Decomposition Time: 12436.2 ms
----------------------------------------
Matrix Size: 768
Serial Matrix Decomposition Time: 31599.1 ms
Vectorized Matrix Decomposition Time: 22726.7 ms
Multi-Threaded Matrix Decomposition Time: 15303.8 ms
----------------------------------------
Matrix Size: 832
Serial Matrix Decomposition Time: 40160 ms
Vectorized Matrix Decomposition Time: 30389.5 ms
Multi-Threaded Matrix Decomposition Time: 19021.6 ms
----------------------------------------
Matrix Size: 896
Serial Matrix Decomposition Time: 51145.3 ms
Vectorized Matrix Decomposition Time: 35971.2 ms
Multi-Threaded Matrix Decomposition Time: 23198.1 ms
----------------------------------------
Matrix Size: 960
Serial Matrix Decomposition Time: 59323.7 ms
Vectorized Matrix Decomposition Time: 46182.9 ms
Multi-Threaded Matrix Decomposition Time: 27722.1 ms
----------------------------------------
Matrix Size: 1024
Serial Matrix Decomposition Time: 75577.6 ms
Vectorized Matrix Decomposition Time: 54494.3 ms
Multi-Threaded Matrix Decomposition Time: 31570.8 ms
----------------------------------------

*/
